#include <algorithm>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <list>
#include <map>
#include <set>
#include <sstream>
#include <string>
#include <queue>
#include <vector>
#include <cmath>
#include <boost/program_options.hpp>

namespace po = boost::program_options;
using namespace std;

int main(int argc, char * argv[])
{
  std::string ofile;
  double kk, posi, width, bin, posi1, posi2;
  double x0, x1;
  
  po::options_description desc ("Allow options");
  desc.add_options()
      ("help,h", "print this message")
      ("strength,k", po::value<double > (&kk)->default_value (4000.), "the k value")
      ("posi,p",  po::value<double > (&posi)->default_value (0.2),   "the first zero point")
      ("width,w",  po::value<double > (&width)->default_value (0.05),   "the second zero point")
      ("x0,l",  po::value<double > (&x0)->default_value (0.0),   "the third zero point")
      ("x1,u",  po::value<double > (&x1)->default_value (1.2),   "the third zero point")
      ("bin,r",  po::value<double > (&bin)->default_value (.001),   "the bin size")
      ("output,o", po::value<std::string > (&ofile)->default_value ("table_b0.xvg"), "the output table");
  
      
  po::variables_map vm;
  po::store(po::parse_command_line(argc, argv, desc), vm);
  po::notify (vm);
  if (vm.count("help")){
    std::cout << desc<< "\n";
    return 0;
  }
  
  unsigned ndata = (x1 - x0 + 0.01 * bin) / bin;
  bin = (x1 - x0) / double(ndata);
  ndata ++;
  cout << "# print table with bin size: " << bin << endl;
  cout << "# ndata: " << ndata << endl;
  posi1 = posi - 0.5 * width;
  posi2 = posi + 0.5 * width;

  FILE * fout = fopen (ofile.c_str(), "w");
  if (fout == NULL){
    cout << "cannot open file " << ofile<< endl;
    exit (1);
  }
  fprintf (fout, "# table of bond, kk %f, posi1: %f, posi2: %f\n",
	   kk, posi1, posi2);
  for (unsigned ii = 0; ii < ndata; ++ii){
    double xx = x0 + ii * bin;
    double vv = kk * (xx - posi1) * (xx - posi2);
    double ff = - (kk * (xx - posi1)+
		   kk * (xx - posi2));
    fprintf (fout, "%f %e %e\n", xx, vv, ff);
  }
  
  fclose (fout);
  
  return 0;
}
